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Abstract 

We propose a method and algorithm for computing the weighted Moore- 
Penrose inverse of one-variable rational matrices. Continuing this idea, we 
develop an algorithm for computing the weighted Moore-Penrose inverse of 
one- variable polynomial matrix. These methods and algorithms are gener- 
alizations of the method or computing the weighted Moore-Penrose inverse 
for constant matrices, originated in [2H], and the partitioning method for 
computing the Moore-Penrose inverse of rational and polynomial matrices 
introduced in [23] . Algorithms are implemented in the symbolic compu- 
tational package MATHEMATICA. 

AMS Subj. Class.: 15A09, 68Q40. 

Key words: Weighted Moore-Penrose inverse; Rational and polynomial 
matrices. 

1 Introduction 

Let C be the set of complex numbers, C mx " be the set of to x n complex 
matrices, and C™ xn = {X e C mx " : rank(A) =r}. For any matrix A e C mxn 
and positive definite matrices M and TV of the orders m and n respectively, 
consider the following equations in X, where * denotes conjugate and transpose: 

(1) AXA = A (2) XAX = X 

(3M) (MAX)* = MAX [AN) (NXA)*=NXA. 

The matrix X satisfying (1), (2), (3M) and (4N) is called the weighted Moore- 
Penrose inverse of A, and it is denoted by X = A^ MN . Especially, in the case 
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M = I m and N — I n , the matrix X — A\ IN comes to the Moore-Penrose inverse 
of A, and it is denoted by X = A^ . 

For any matrix A £ Q nxn the Drazin inverse of A is the unique matrix, de- 
noted by A D , and satisfying the matrix equation (2) and the following equation 
inl: 

(l fe ) A k XA=A k , (5) AX = XA. 

As usual, C[s] (resp. C(s)) denotes the polynomials (resp. rational func- 
tions) with complex coefficients in the indeterminate s. The m x n matrices 
with elements in C[s] (resp. C(s)) are denoted by C[s] mx " (resp C(s) mx "). 
By I is denoted an appropriate identity matrix. 

We observed three different directions in the symbolic computation of gen- 
eralized inverses: 

A) extensions of Leverrier-Faddeev algorithm, 

B) methods based on the interpolation, and 

C) methods based on the Grevile's recursive algorithm. 

A) Computation of the Moore-Penrose inverse of one variable polynomial 
and/or rational matrices, based on the Leverrier-Faddeev algorithm, is investi- 
gated in [TJ 13 LHJ [TJ1 LH3 [5S] . These papers are based on the paper [5] . Imple- 
mentation of the algorithm from [12] in the symbolic computational language 
MAPLE, is described in 11 . An algorithm for computing the Moore-Penrose 
inverse of two- variable rational and polynomial matrix is introduced in |16) . A 
quicker and less memory-expensive effective algorithm for computing the Moore- 
Penrose inverse of one- variable and two- variable polynomial matrix, with respect 
to those introduced in [TJ] and [TB], is presented in [TJ]. This algorithm is effi- 
cient when elements of the input matrix are polynomials with only few nonzero 
addends. 

Continuing the algorithm of the Leverrier-Faddeev type for computing the 
Drazin inverse of constant matrices, established in [6], a representation and 
corresponding algorithm for computing the Drazin inverse of a nonregular poly- 
nomial matrix of an arbitrary degree is introduced in |10j . [221 125] . Bu and 
Wei in [3| proposed a finite algorithm for symbolic computation of the Drazin 
inverse of two- variable rational and polynomial matrices. Also, a more effec- 
tive three-dimensional version of these algorithms is presented in the paper [3] . 
Implementation of this algorithm in the programming language MATLAB is pre- 
sented in [3]. 

A general for of the Leverrier-Faddev type algorithms is introduced in [24] . 
This algorithm generates the class of outer inverses of a rational or polynomial 
matrix. 

B) An interpolation algorithm for computing the Moore-Penrose inverse of a 
given one- variable polynomial matrix, based on the Leverrier-Faddeev method, 
is presented in [20 . Algorithms for computing the Moore-Penrose and the 
Drazin inverse of one-variable polynomial matrices based on the evaluation- 
interpolation technique and the Fast Fourier transform are introduced in [15] . 
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Corresponding algorithms for two-variable polynomial matrices are introduced 
in |27j . These algorithms are efficient when the input matrix is dense. 

C) Grevile's partitioning method for numerical computation of generalized 
inverses is introduced in [7]. Two different proofs for Greville's method were 
presented in [4], [29]. A simple derivation of the Grevile's result has been given 
by Udwadia and Kalaba [26] . In [8] Fan and Kalaba used the approach of deter- 
mination of the Moore-Penrose inverse of matrices using dynamic programming 
and Belman's principle of optimality Wang in [28] generalizes Grevile's method 
to the weighted Moore-Penrose inverse. Also, the results in [28] are proved using 
a new technique. 

In [21] the Greville's algorithm is estimated as the method which needs more 
operations and consequently it accumulates more rounding errors. Moreover, 
it is well-known that the Moore-Penrose inverse is not necessarily a continu- 
ous function of the elements of the matrix. The existence of this discontinuity 
present further problems in the pseudoinverse computation. It is therefore clear 
that cumulative round off errors should be totally eliminated. During the sym- 
bolic implementation, variables are stored in the "exact" form or can be left 
"unassigned" (without numerical values), resulting in no loss of accuracy dur- 
ing the calculation [12) . 

An algorithm for computing the Moore-Penrose inverse of one- variable poly- 
nomial and/or rational matrices, based on the Grevile's partitioning algorithm, 
was introduced in [23 . An extension of results from [33] to the set of two- 
variable rational and polynomial matrices is introduced in the paper [19] . 

In the present paper we extend Wang's partition method from [28) to the 
set of one-variable rational and polynomial matrices. In this way, we obtain 
an algorithm for computing the weighted Moore-Penrose inverse of one- variable 
rational and polynomial matrices. The paper is a generalization of the paper 
[28] and a continuation of the paper [23] . 

The structure of the paper is as follows. In the second section we extend 
the algorithm for computing the weighted Moore-Penrose from [28] to the set of 
one- variable rational matrices. In Section 3 we give the main theorem and adapt 
this algorithm to the set of polynomial matrices. Several symbolic examples are 
arranged in fourth section. In partial case M — I mi N — I n we obtain the usual 
Moore-Penrose inverse, and then use test examples from 32 . In the last section 
we describe main implementation details. 

2 Weighted Moore-Penrose inverse for rational 
matrices 

Greville in [7] proposed the partitioning algorithm which relates the Moore- 
Penrose pseudoinverse of a constant matrix R augmented by a vector r of ap- 
propriate dimensions with the pseudoinverse RJ of R. Wang and Chen in [28] 
generalize Greville's partitioning method. They obtained an algorithm for com- 
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puting the weighted Moore-Penrose inverse, and give a new technique for its 
proof. This method is also suitable for the weighted least-squares problem. 

By Ai(s) we denote the submatrix of A(s) G C(s) mx ™ consisting of its first 
i columns: 



Ai(s) = Aj_i(s) | di(s) , i = 2,...,n, M{s) = d\(s) 



(2.1) 



where cii(s) is the i-th column of A. 

In the sequel we consider positive definite matrices M(s) G C(s) mxm and 
N(s) G C(s)" x ™. The leading principal submatrix Ni(s) G C(s) tx ' 1 of N(s) is 
partitioned as 



Ni{s) 



iVi-i(a) 
2*( s ) 



2,...,n. 



(2.2) 



In the following lemma we generalize the representation of the weighted Moore- 
Penrose inverse from [28] to the set of one- variable rational matrices. 

For the sake of simplicity, by Xi(s) we denote the weighted Moore-Penrose 
inverse corresponding to submatrices Ai(s) and iVj(s): Xj(s) = Ai(s)\ IN . , for 
each i = 1, . . . , n. 

Lemma 2.1. Let A(s) G C(s) roxn J assume that M(s) G C(s) mxm and N(s) G 
C(s) nx " are positive definite matrices, and let Ai(s) be the submatrix of A{s) 
consisting of its first i columns. Assume that the leading principal submatrix of 
N(s), denoted by Ni(s) G C(s) iX % is partitioned as in (jg.gp . 
In the case i = 1 we have 



X 1 (s) =ax(sV = { K( s ) M ( s ) a i( s )) l at(s)M(s), Oi(«) ^ 0, 
\ aj(s), oi (s) = 0. 

for eac/i « = 2, . . . , n Xi(s) is equal to 
x (s) J (*(*) + (I- X<-i (s)li-i(s)) iVr\( s )^( s ))6|( s ) 

where the vectors di(s), Ci(s) and b*(s) are defined by 
di(s) = Xi_i(s)di(s) 

c,(s) = a,i(s) - Aj_i(s)d»(s) = (/ - Aj_i(s)Xj_i(s) ) Oj(s). 



(2.3) 

(2.4) 

(2.5) 
(2.6) 



(2.7) 



(c*(s)M(s)c t (s)y 1 c*(s)M(s), a(8) ? 

K\ 3 ) (d?(«)JVi-i(a) - Xi-i(s), Ci(a) - 0, 

and where in <\2. 7] ) is 

o"i(s) = n«(a) + d|(s)iV i _i(s)d i (s) - (d*(s)h(s) + l*(s)d l (s)) 

-l*(s) (/ - iVr\ ( S )^( s ). (2.8) 
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Proof. The proof is the same as for constant matrices, presented in [55]. □ 



The following lemma is a simple extension of the well-known result in the 
literature. 

Lemma 2.2. Let A(s) be a partitioned matrix which is nonsingular, and let the 
submatrix An(s) also be nonsingular. Then 



where 



A(s) 



B n (s) 
B 12 (s) 
B 21 (s) 
B 2 2{s) 



An(s) A 12 (s) 
A 2 i(s) A 22 (s) 



Bn(s) B 12 (s) 
B 21 (s) B 22 (s) 



An(s)- 1 + A n (s)- 1 A 12 (s)B 22 (s)A 21 (s)Ai 1 (s) 

-A n {s)- l A 12 {s)B 22 (s) 

-B 22 (s)A 21 (s)A 11 (sy 1 

(A 22 (s) - A^Auis^Auis))- 1 . 



The following lemma is a generalization of known result from 
of rational matrices. 

Lemma 2.3. Let Ni(s) be the partitioned matrix defined in 
Ni(s) and Ni-i(s) are both nonsingular. Then 



Wi-i(s) k(s) 

l*(s) Uji(s) 



Ei-i(s) Ms) 
ft(s) 9u(s) 



where 



(2.9) 

- 1 (2.10) 
(2.11) 
(2.12) 
(2.13) 

to the set 
ssume that 

(2.14) 

(2.15) 
(2.16) 
(2.17) 

= Ni-i{a), 



In view of Lemma 12.11 we present the following algorithm for computing the 
weighted Moore-Penrose inverse of a given one- variable rational matrix. 

Algorithm 2.1. Input: rational matrix A(s) £E C(s) mxn and positive definite 
matrices M(s) G C(s) mxm and N(s) £ C(s) nxn . 

Step 1. Initial value: Compute X\{s) = a\{sy defined in §2.3\) . 

Step 2. Recursive step: For each i = 2, . . . , n compute Xi(s) performing 
the following four steps: 

Step 2.1. Compute d,(s) using (|S.5jl . 
Step 2.2. Compute ct(s) using (|g.6]) . 
Step 2.3. Compute b*(s) by means of §2J\ and (TO|) . 
Step 2.4- Applying \2.4\ compute Xi(s). 
Step 3. The stopping criterion: i — n. Return X n (s). 



9ii{s) 

Ms) 

E l - l {s) 



(n^-ltmr^sMs))- 1 
-gu(s)Nr_\( s )h(s) 

N-M^ + g^i^M^f*^)- 



Proof. The proof immediately follows from the substitutions An(s) - 
A 12 (s) = k{s), A 21 (s) = li(s)* and A 22 {s) = na(s) in Lema\E^ □ 



6 



M.B. Tasic, P.S. Stanimirovic, M.D. Petkovic 



Let iVj(s), defined in (|2.2p . be the leading principal submatrix of positive 
definite matrix N(s) £ C(s) nxn . The following algorithm, based on Lemma 
12. 31 computes the inverse matrix (s) £ C(s) ixi . 

Algorithm 2.2. Compute N~ 1 (s). 

Step 1. Initial values: N± (s) = n^(s) 

Step 2. Recursive step: For % = 2, . . . , n perform the following steps: 
Step 2.1. Compute gu(s) using <\2. 15$ . 
Step 2.2. Compute f%{s) applying \2.16\ . 
Step 2.3. Compute Ei-i(s) according to (\2. 17$ . 
Step 2.4. Compute N j ^ 1 (s) using ( [gjg ). 

Step 3. Stopping criterion: for i — n the output is the inverse matrix 
N- 1 {s) = N-\s). 

3 Weighted Moore-Penrose inverse for polyno- 
mial matrices 

Consider the matrix A(s) £ C[s] mxn given in the polynomial form with respect 
to unknown s: 

A(s) = A-l + A 2 s + • • • + Aqsi- 1 + A q+lS i = ^ A l+1 s l (3.1) 

j=Q 

where Ai, i = 1, . . . , q + 1 arc constant m x n matrices. 

Theorem 3.1. Consider an arbitrary polynomial matrix A(s) £ C[s] mx ™ given 
by \3. 1\ and the following polynomial forms of positive definite matrices M{s) £ 
C[s] mxm and N(s) £ C[s] nxn : 

m g n q 

M(s) = M i+ is\ N(s) = J Ni+is*. (3.2) 

i=0 i=0 

Transcribe i-th column of A(s) by 

q 

ai(s) = y^ y a it j + is j , l<i<n, (3.3) 
3=0 

where djj+i, < j < q are constant m x 1 vectors. Also, denote first i columns 
ofA(s) by 

A i (s)=J2Aj+i sj , l<*<n, (3.4) 
3=0 

where Ajj+i, < j < q are constant m x i matrices. 
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In the partition <\2. e 2ty of the leading principal submatrix Ni{s) £ C[s] ixi of 
N(s), we use the following polynomial representations: 

n q 

n q n q E Ni-ij + iS j 

n li {s) = '^ri l . : j + is j , k(s)=y^^Li t j + is j , Njli(s) = ^ • (3.5) 

E ^-u+i^ 

TTien the following algorithm computes the weighted Moore-Penrose inverse 

Algorithm 3.1. 

Step 1. Initial values: 

Compute Zi t j+i, 0<.j<q\=q + m q and Y\_j + \, < j < p% = 2q + m 9 as in 

{E a i j-fe+i-Mfe+i, < j < <7 + to 9 , oi(s)^0, 
(3.6) 
a ij+i=°> a x (s)=0. 

E E aij-jfe-r+i-Wk+iai.r+i, 0<j<2g+m g , ai(s)^0, 
r=0fc=0 (3. 7 ) 

1, ai (s)=0. 

S^ep £. Recursive step: 

For 2 <i <n perform Step 2.1 , Step 2.2, Step 2.3 and Step 2.4. 
Step 2.1. Compute di.j+i: < j < q + by means of 

j 

<k,j+l = E Zi-l.j-k+iai.k+l, < j < qi-i + q. (3.8) 
fe=0 

Step 2.2. Compute Ci,j+i> < j < q\-\ + q using 



3 

c i,3+i=^2( a i,3-k+iYi-i,k+i - j4j-ij-fc+i4,fc+i), 0<j<gj_i+g. (3.9) 

fc=0 

where 



q,-! = max{p l _ i,<?i-i (3.10) 
S'iep 2.3. If Cij+i t^O /or some j, compute Vij+i and Wij+i by means of 
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3 3-r 

= ^ ^ Yj-l,j-k-r+lC* k+1 M r+ i, 
r=0 k=0 

0<j<bi = Qi-i + q + Pi-i + m q , 

3 j-r 

^ ^ clj- k -r+l M k+lCi,r+l, 



(3.11) 



W, 



i,j+i 



(3.12) 

r=0 fc=0 

0<j <li = 2qi-i+2q + m q . 
In the case c ij+1 = for each j, compute Vij+i and Wij + i in this way: 



j j-t j-t-r 

Vj,j+i = y^ yi y^ Ajj_fc_ r _ t+ idj it+ i./v ; ;_i ir+ i.Zj_i i fc + i 

t=0 r=0 fc=0 
— ^i,j-k-r-t+lL iit+ iYi-\ t r + iZi-i i k+l, 

0<j<bi = 2pi-i + n q + qi-i + &_i + n„ 

j j-r 

+ 1 — y^ ^ Aj,j-k-r+lYi-l,k+lYj-l,r+l, 



(3.13) 
(3.14) 



r=0 k=0 

0< j < &i = 2gi_i + n g + max{n, + n q ,n q } + 2pj_i, 



where 



A i,j+l = yi-l,j-fc-r+l^i-l,fc+l^'i-l,r+ia J ' ! < j < <5, = 2pi_i + n 9 (3.15) 

r = fc = 
J j-t j-r— t 

Aij + l = ^ (nij-fc-r-t+l^t-l.fc+l^i-l.r+l + d*j_ fe _ r _ t+1 -/Vi_l i fc + lfifi, r +l 



i=Or=0 fc=0 



lJ -fc-r-t+l-^<,fc+l I i-l,'-+l ~ i* j j_ fc _ r _i + idi,fc + l^-l,r+l)A r i-l,t+l (3.16) 



- i *,j-fc-r+l¥'i,fc+l i »-l,r+li < j < 5 q = 2gj_i + ?iq + max{n 9 + n q ,n q } 



Then compute 



Zi,j+i, 0<j<q t , Y ifj+1 , 0<j<p l 



as it is defined in 



Zi.-i 



3+1 



< j < ft, i > 2 (3.17) 
= £^j-fc+iWi,fc+i. 0<j< Pj , i>2, (3.18) 



J2 i>i,j-k+lVi,k+l 
k=0 



k=0 
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where in ( fXT7| ) and (fX7^|) 



J 3—r 

= J]] ■Z'i-ljj-fc— r+l^V»-l,fc+l Wi,r+1 
r=0 fc=0 

— rfij-fe-r+l-^i-l,fe+l^i,r+l — ViJ-fc+l^.fe+li (3.19) 

0<i<g l; 

i j'-t 3-t-r 

=y (^-i,2-fc-r-t+x^-i,fc+i _ ^-i,j-fc-»— t+i-^x-i,fc+i) 

i=Q r=0 fe=0 

xJVi-i^iA-.jk+i, (3.20) 
< 3 < 3t-i + + »V 

i 

V'irf'+i = y^^-i.j-fe+i^-i.fc+i; < j < pi-! + n q (3-21) 



fc=0 



and 



qi = Qi-i +q + max.{n q + nq,n q } + max{bi,bi}, q 1 =q + m q (3.22) 
Pi = Pi-i +n q + bi, pt = 2q + m q . (3.23) 



Step 2.4. Compute 

qi 

E Z i,j+l sJ 

X *( s ) = TT . ( 3 - 24 ) 

E 

where qi and pi are defined in \3. 22\ and \3. 23$ , respectively. 

Step 3. The stopping criterion is i — n. In this case the result is the weighted 
Moore-Penrose inverse A(s)^ M N = X n (s). 

Proof. If ai (s) = 0, in view of the second case in (|2.3p we have 



*i00 = = 5>i J+1) y 1 ( S ) = i. 



If ai(s) 7^ 0, in accordance with the first case in (|2.3[) we have 
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Xi(s) = (aJ(s)M( S )a 1 ( S ))- 1 at(s)M( S ) 



-l 

9 



Y a *,j+i sj Y M 3+^ sl Y a hj+i sj Y a h+i s3 Y M j+i s ° 

J=0 j=0 3=0 J j=0 j=0 

'2q+m q j j-k \ 1 q+rn q j 

Y Y^^J-k-r+lMk+iai.r+lS 3 \ Y ^2 a l,j-k+l M k+lS J 
j=Q r=0 k=0 J j=0 k=0 

q+m q j 

E E 

3=0 k=0 



E E<i- fe+ iM fc+1 ^ 



2g+m, j j-r 

E E E a* j -k-r+l M k+l a l.r+lS j 
j=0 r=0 k=0 



Therefore, X\(s) is the partial case i = 1 of (|3.24p . where the matrices Zij+i 
and Yij + i are defined in (|3.6[) and (|3.7p , respectively. 

For each i = 2, . . . , n it is reasonable to calculate matrices -Xj(s) in the form 
(|3.24p . for appropriate matrices Zij+x, Yi,j+i and appropriate upper bounds qi 
and pi . 

Direct calculation in (|2.5|) . i.e. 5£ep U.J of Algorithm 2.1 yields the following 



E ^-lj + lS 3 q 

di(s) = Xi-x(s)a,i(s) = |^ y^aj : fc +:L s fc 

E fc =° 

E ( E Zi-ij-k+iai_k+i)s 3 

j=0 k=0 

~ Pi-i " 

E *i-i,j+is j 

i=o 

Then di(s) can be represented in the form 

9,-1+9 

E d i,j+i s3 

d ^ s ) = pT~r . ( 3 - 25 ) 

E Yi-WS j 

3=0 

where the matrices djj+i are defined by (|3.8p . 

Consider ()2.6)) . i.e. Step 2.3 of Algorithm 2.1. Since the first j — 1 columns 
of A(s) can be represented in the polynomial form 

9 

Ai-i(s) = Y Aj-i.j+is 1 

3=0 
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for appropriate m X (t - 1) constant matrices Ai_ij-|_i(s), in view of Q3.3P and 
(|3.25|) we obtain 



Ci(s) = ai(s)-Ai-i(s)di(s) 



qi-i+q 

E 

J=0 



E [ E a i,i-fe+i^-i,fc+i Is 3 — (E Ai-i,j-k+idi,k+i) 

j=0 \k=0 J j=0 k=0 

Pi-l 

E 

j=0 



Finding a maximum between the upper bounds q + pi-\ and 2q + qi-\ in 
the last identity, we have 



94-1+9 j _ 

E E ( a i,j-k+iYi-i.k+i — Ai^ij-k+idi^+^s 3 

I n 3=0 k=0 

E Yi-i,j+is j 



where q~i-\ is defined in (|3.10p and shorter polynomial matrix is filled by appro- 
priate zero matrices. 

Therefore, Cj(s) can be represented in the form 



qi-i+q 

E Cij+is 3 

I \ 3=0 
Ci(S) 



Pi-l ■ 
E Yi- 1J+1 si 

3=0 



where Cj are matrices of the form (I3.9p . for each 0< j + 9. 

Observe now Step 2.3. of Algorithm 2.1, i.e (JHU). 
If Cjj+i 7^ for some j, then Ci(s) ^ and b*(s) is equal to 
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b*(s) =(c*(s)M(s)c i (s))- 1 c*(s)M(s) 



9.-1+9 



E 



Si-i+g 

E Ci, i+ iS j 

£^Ws_? 

E ^-u+i^' i=° E ^-u+i^ 



Pi-i 



Pi-i 



-l 



9.-1+9 

E 

i=0 



E C ij'+l S "' rn q 



Pi-l 



E Yr-w^- E E ct ._ fe+1 M fc+lS ^ 

■>=n j=0 k=0 



3=0 



2q i - 1 +2q+m q j j-r 

E EE c* 1 _ k _ r+1 M k+1 c l , r+1 s3 

— n » — n u — n 



j=0 r=0 fc=0 

9i-l+9+Pi-l+m, j j— r 



+Pi-l+m, j j-r bi 
E EE ^-lj-fc-r+lC* fc+1 M r + iS 3 E ^J + 1 S 

3 =0 r=0 k=0 j=0 



2q i - 1 +2q+m q j j-r 

E EE cf •_ k _ r+1 M fc+ iCi,r+i* i v Wi 

3=0 r=Ok=0 l^n 



sJ 



where Vij + i and Wij+i satisfy (13. lip and (|3.12j) . respectively. 

If Ci,j+i = for all j, then Cj(s) = and b*(s) is defined in the second case of 
(|2.7|) and in (|2.8j) . In order to compute (S i " 1 (s), we firstly generate the following 
intermediate value, which will be used later: 



CT4 (s) = (I - X i ^ 1 {s)A i ^ 1 (s)) Nr_\(s)li(s) 



( n-i 

'-^ E^-w' 

\ 3=0 



= E L i.J + l s3 

3=0 



/ E ffi-i, j+ iai 

/ 3 = 



9 z — 13 , 1g+«5 j 

E E Yi-l,j-k+lli-l.k + l — Zi-l,j-k + lAi-i,k + lS J ^2 E ^ i-i-,j-h + lLik + lS J 
3=0 k = 3 = k=0 



Pi-1 

E 

3 = 



.S'J 



3=0 



+ + 3 j — tj — t — r _ 

E E E E (Yi — 1,3'— fe — ' — t+ll%— l,fe+l — 1,3 — fe — 1 — t+lAi — x,k+l)Ni — l,r+l-^i,fe+l 

3 = t = r=0 fc=0 

Pi-1+5?, j _ 

E E 3-^ + 1-^4-1, fc+is3 

3=0 fc=o 

"!i-l+n : 5+"Qr 

E Vi^ + i* 3 

3 = 

Pi-l+Ttq 

E ^i,3+l sJ ' 



In the last identity ifij+i and tpij+i and <5i are defined by (13.20P and Q3.2ip . 
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Now, Si(s) is equal to 

6i(s) = riii(a) + d*(s)N i - 1 (s)d i (s) - (d*(s)/,(s) + Z*(s)d s (s)) - /*(s)cr,(s) 

n q E + n, E diJ + lS 3 



3=0 3=0 
Qi_l+Q g*_l+g Qi-l+^q+^g 

E ^.J + l^ », E d .J + l sI n, E fi,0 + lS 3 

Pi-1 Z ' 3 + 1 ^ + 1 Pi-1 ^ + 1 p ,+= 

3 = 3 = — „ 

2q i _ 1 +n„ j j_ 



sr V 3 = 3 = ' W v 3 = Pi-l+ng 

j?0 + 1 3?0 Yl - 1 ' J + 1 jg, ^,3 + 1^ 

-1+n, j J-r 

E E E "i,j-i-r+l^i-l.* + lli-l,r + l + ,-fc-r+l N i- 1 ,k + 1 di , r +l S J 

3 = r=0 fc = 

E -1 E yi-i,i-*+iVi_i,* + isJ 

3=0 fc = 

— l+^g 3 3 — r 

E E E ^ i ,-t- r+ li.,Hl i 'i-l,<'+l + I *j-t-,+ l''i,Hl ! 'i-l,r+l s3 

3=0 t-=0 fc = 



2 Pt-l 3 

E E y;-i,3-fc+i^-i,fc+i^ 

3=0 fc=0 

E E l* ] _ k+1 ipi, k+ is 3 

3=0 fc=0 

E V-ij+i* 3 ' 

3 = 

2q i _ 1 +n q + n q j j-tj-n 1 _ 

E E E E (fii,j-k-r-t+iY i - lik+1 Yi- 1:r+1 + d* ■ k _ r _ t+1 N i _ 1 . k+1 di tr+1 )N i ^ ltt+1 s : > 

3 = t=0 r=0 fc = 

2p i _i+W, 1 3 3-r = 

E EE ^i-l.J-fc-r+lii-l.fc+lJVi-l.r+lS^ 
3 = r = 0fc = 

2q i _l+nq+—q j j-t j-r-t _ 

E E E E (d* j _ k _ r _ t+1 L i , k+1 Y i - 1 , r+1 + L* j _ k _ r _ t+1 d iik + 1 Yi- lir+1 )N i - ltt+1 s : > 

3=0 t=0r=0 k=0 

2pi-l+rTq j j — r _ 

E EE y i -i,j-*-r+iyi-i,fc+iJv i _i, r+ isJ 

3 = r=0 fc=0 

9i_l+2n, , 3 j_ r 

V E V L* ID- l.,K i j_,S 3 ' 

3 = r=0 fc = 



2p i _ 1 +n q j j_ r _ 

E EE y»-l,3-fe-r+l'5 / i-l,fe+l-'V»-l,r+lS ;i 
3 = r=0 fc = 



5q j j — t j — r — t — 

E E E E ("ij-fc-r-t+iy-i.fc + iy-i.r+i + d* ,_ fc _ r _ t + 1 -'V i _ 1 j c + 1 d iir+ i)Ar i _i >t+ is 3 

3 = 0t=0r-=0 fc = 

2 Pz-l+~g j 3-r 

V 2 E Yi-i, j -fc-r+iV i _i,*+iAr 4 _i, r+ isi 

j = r=0 fc = 



5 Q j j — t j — r — t — 

E E E E (dr,j_fc_ r _t + i^i,fc+iVi-i,r+i +^* J _ fc _ r _ £+ irfi,fc+i"s / i-i,r+i)Ari_i )t+ is J 

j=0t=0r=fl fc=0 

2 Pi-l+~g j j-r — 

E EE ^"i-l,J-fc-r+l^i-l J fc + lJV'i-l 1 r+ia J 
j = r=0 fc = o 

5 Q j j — t g Q — 

E E E ^L-fc-r+i^.fc + lVi-l.r+l^ E A>i,j + lS 3 

j=Qr=0k=0 ' j=0 



P 'e q E E^i-u-fc-r+i^-Lfc+i^-i^+i^ E A,, j+ is^ 

j = r=0 fc = j = 
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Therefore 



Si(s) 



-i = ^ 



where A, j+i and Aj J+ i are defined in (|3.15[) and (|3.16l) . respectively. 
Now, in accordance with the second case of (|2.7|) . &i(s) is equal to 



b*(s) = S-^s) (dUsJN^is) - l*(s)) 



52 A iij+ isJ 
j=o 

?, = 

2 A ijj+ isJ 

3=0 



E ^i.j + l sJ n q n q 

— E N i-u+i sj - E 

3=0 



. ^ ^-1,3 + 1*"=° 
\ 3=0 



5] Zi- i,3+i * J 

3=0 

Vi-1 
E Vi-l,3 + lS j 
/ 3 = 



Oq / 3 

E Ai,j+is J '. E E ^-fe+i^i-i.fc+i - L lj-k+i 

3=0 3=0 \k=0 



U-l 

3=0 



Sg+2pj_i / 3 j-r_ \ 

E (EE ^i,j-fc-r+l^»-l,fc+l^»-l,f+l I 
3=0 \r=0fc=0 / 

*<j+9i-l+<3s-l+™<} j j — tj — t — r 

E E E E ^■ij-k-r-t+ldf k+1 Ni-l !r+ lZi-lj + lS 3 

3=0 t=0r=0 fc = 

5 ? +2pi_i / j 3-r_ \ 

E (EE Aij-fc-r + l^-l.fc + l^-l.r+l I * J 
j=0 \r=0fc=0 / 

E E E E ^i.j-k-r-t + lL* k + 1 Y i - ltr + 1 Zi- 1 . t + 1 S^ 

j=0 i=0r=0 fc=0 

5,+2p i _ 1 / 3 3-r_ 

E EE ^i,j-*-r+l^i-l,fc+l5»-l,r+l I « J 

3 = \r=0fc = 

E ^.3 + lsJ 
3=0 



E Wi,3+1^ 
3=0 



It is not difficult to verify that in the last expression Vij+i and Wi j+i satisfy 
(|3~T3f and (|3TT4"|) . respectively. 



Finally, using (|2.4j) of Algorithm 2.1, we obtain 
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Xi{s) 



Xi_i(s) - (di(s) +cn(s))b*(s) 



9i-l 

E Zi_ 1J+1 s^ 

3 = 

Pi-1 

e n-i,i+i» j 



E d i,3 + l sJ £ fi.J + l^ 



Pi — 1 

e n-i,i+i« 

5=0 



+ 



Pi-l+^g 

E V'i.j+isJ 

3=0 



/ E w*,i+i^ 

/ 3 = 



3=0 



E w s . 3+lS 3 

3 = 



9i-l 

E Zi-lj + is 

3=0 

Pi-1 

E Y *-1 >3 '+1< 



?i — j _ Qi-l + ^g+^g 6^ 

E E d i j _ k+1 N i _ l k+1 .i+ E *>i,3+i si E 

3=0 fc = 3=0 3=0 



P»_l+»g 

E ^i,fc+i sJ 



E W i.j + l s 
3 = 



3—0 



E 

3=0 



Qi-l 

E Zi-ij+i* 

3=0 

Pi-1 

E Vi-1,3+1* 3 
3 = 



1i— \+1+ n q j <Zi — l+'^q^ 

E E i, lI -Hi» i -u + i» ] + E 
3=0 fc=o ' T 3=0 



J = 



E „ e 4,i + l' 



E E ^j-fc+i w i,Hi s 

j = fc=0 



3=0 



J=0 



E e i,3 + l S 
3=0 



Pi — l+™g+&i j 

E E tfi j-k+l w i,k+-L e 

j = fc=o 

Pi-l+^q+^i j 

E E *i,,'_fe+l Vi^fc+l" 

3 = fc = 

Pi_ 3 

E E ♦ij-Hi^.H!" 

3=0 fc=o 



Pi — l + n q + b i 3 

E E +, J -Hl»'i,i + l J 

3=0 fc=o 

E V i,3 + 1 sJ 
J=0 



J = 



E °i,3 + l s 
3 = 



Pi — l+n q + b i j 



j=0 k=0 



E V i,3+1 6 
j=Q 



3 = 



Pi-l + "9 

E ^i,3 + l S 
3 = 

Pi-l + ™9 

E l*ij+l« 
3 = 



E n z i,3+1 = 

Pi 

E n,3+ia 
3 = 



where is defined in (|3.19p . 

Finally, we obtain the polynomial representations for Zjj+i and lij+i as in 

dsmi-dsisi) 
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In accordance with Lemma 2.1, the weighted Moore-Penrose inverse for given 
t 

M,N 



matrix is j4(s)j^ n = X n (s), which completes the proof. □ 



The next algorithm is a generalization of Algorithm \2.2\ and computes the 
inverse matrix N~ 1 (s) in a polynomial form. 

Theorem 3.2. Let the leading principal submatrix Ni(s) of the positive definite 
matrix N(s) is partitioned as in (12. 21 , and assume that nu(s),li(s),N~\(s) 
possesses the polynomial representation (|ff.5[) . Then the following algorithm 
computes the inverse matrix N~ 1 (s). 

Algorithm 3.2. Input: positive definite matrix N(s). 
Step 1. Initial values: 

Ni, j+ x = 1, = n hj+ i, 0<j< n q . (3.26) 

Step 2. Recursive step: For 2 < i < n perform Step 2.1-Step 2.4.: 
Step 2.1. Compute 

Gij+i = < j < g q = % (3.27) 

j 

P»,i+i = X! ™i,j-k+iNi-i,k+i, < j < n q + n q 

k=0 
3 3-r 

= ^2^2 L i,j-k-r+l N i-l,k+l L i,r+l, < j < 2n q + n q 

r=0 k=0 

Gij + i = Vi,3+\ - < j <g q = 2n q + n q (3.28) 

where Pij+i is padded by zeros from n q + n q up to upper bound 2n q + n q . 
Step 2.2. Compute 

3 

Fi, J+ i = - } j A r i-ij-fe + iLi, fc+ i, < j < f q = n q + n q (3.29) 

fc=0 

%, J+ i = Sij+i, < j < J q = % (3.30) 
Step 2.3. Compute 

3 3-r 

Ei-l.j + l = / ] Nj-ij-k-r+1 G't.fe+l-P'i.r+t 

r=0 k=0 

+ Ni-ij-k-r+lFi t k+lF i _ lr+1 , (3.31) 

<j < e q = ma,x(n q + g q + f q ,n q + 2f q ), 
Ei-ij+i = Ni-ij-k-r+iGi^k+iFi_i }r+ i, (3.32) 
<j <t q = % + g q +7 q - 
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Step 2-4- Generate 



j j~r __ __ j j—r^ = __ 

E E Ei,j-k-r+lFi,k+lGi, r +l E E Ei,j-k-r+lFi,k+lGi,r+l 

r=Ofc=0 r=Ofc=0 

j j—r^ = __ j j—r __ __ 

E X ^'i.i-fc-r+l-F'i.fe+lG^r+l E X -E'ij-fc-r+l-F^fc+lG.Lr+l 

r=Ofe=0 r=Ofc=0 



(3.33) 

< j < n q = max§, + 7 g + e g , V q +7 q + V q +7 q + 5g, 9 q +7 q + ^} 

J 3-T 

Ni.j+i = J5^ l j_jfe_ r +xi ? 'i,fe4-iGi,r+l) 



r -0 fe=0 



(3.34) 



< j < n g = g q + f q + e q . 



Step 3. Stopping criterion: for i = n the inverse N (s) = N n 1 (s) is 
to 



E N n , j+1 si 

3=0 

n, _ 

E ^ J+ is J 



Proof. It is not difficult to verify that (|3.26[) follows from 

1 



E riij+isj 

3=0 



(3.35) 



Also, (pOTl) . (|3~^5)l , (15301) . ([5HD and (15321) follows from the following. 

Using (|2. 15[) we have 



9 q 9q _ 

3=0 3 = 

( ?*- \ _1 

I 7\r „i 



E^ni.j+is^ - E^ l »*,j+i" ■= 

3=0 2=0 



3=0 



V 



E Ni- lt j +1 si 

3=0 



E^ Li >j+ isJ 

3=0 



E Ni-ij+is* 

3=0 



"q+n q j _ 2n q + n q j j- r 

E E fti,3-k+l N i-l,k+l s:i - E E E ^.'j-lb-r+l^i-l.Ji+lii.r+lS' 



j=0 fe=0 



j=0 r = 0fc = 
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An application of (|2.16p gives 
7 7 

AM = (E^,,+i^')/(E^, J+ i^) = -9^(s)N-_\(s)h(s) 

3=0 3=0 



- E G i.j+i sJ E Ni-ij+i^ ,i„ 

3=0 3=0 V- j 
= 2^ L i.3 + l S 



2n q +n q _ 

E 

3=0 j = 

In view of (|2.17l) one can verify the following: 



™q+™g j 

' E E ^i-l.j-k + lLi t h + lS 3 
3=0 fc=0 

2n q +n q _ 

E Gi-lj + lS' 
3=0 



Si-iW = (E B ^+i si )/(E E ^+i s3 ) = Jv il\M+s« 1 MAM/fM 

3=0 3=0 

^g 17g _ /g /g 

E Ni-u+iJ E Gij+i^' E ^,3+i* J E Kj+i s3 

3=0 , 3 = 3 = 3=0 



n q _ 9 q f j- 

ZNi-ij+isi E^,3+i^' 

3=0 3 — U 



>,3 + J- a " Z-j i,3+l 4 
3=0 3=0 



max(n q +g q + 2g q ,n q + g q +2g q ) j j_ r ^ _ 

E EE 'Vi-lj-fe-r + lGi^+lFi^+i 

3=0 r=0 fe = 

n q +9 q +f q j j-r _ __ 

E EE ^t-l,i-fc— r+lG^fc+lFijr+l* 7 

3=0 r=0 fc = 

max (rig + g q +2g q ,7T, + g q + 2g q ) _ 

E N i-l,j-k-r + \Fi k + lF iyr + lS 3 
3=0 

= g + 9q+7g 3 j-r— _ — 

E EE N i-lJ-k-r+lGi^+lFi.r+lS- 1 
3=0 r=0fc=0 



Using (|2.2[) we finally get the inverse 



E;_l(s) /i( S ) 

fT(s) gu(s) 




9q-*rfq+ e q 

E 

3=0 


3 3— r 

E E 

r=0 k = 


B s 




-k- 








9q + fq + e<l 

E 

3=0 


3 3 — r 

£ £ 

r=0 k=0 


E, 




-k- 


r+lF 




,r+l 


Sg + /g+^g 


j j-r 














Hq + f q + e q 


j j — r 














E 


£ £ 


~E, 




-k- 


r+1?, 






E 


£ £ 


~E, 




-k- 




,k + lGi 




3=0 


r=0 k = Q 










3=0 


r=0 k=0 










9 q + f q +e q 
















9g + /q + eg 


j j — r 














E 


£ £ 


£ ; 


J 


-k- 


r+l F i,Hl G > 


,T+1 


E 


£ £ 


~E< 


>j 


-k- 




,k + \G % 




3=0 
















3=0 


r = Q k = 












9 q + f q +e q 

E 

3=0 


j 3 ~ r 

£ £ 

T = fc=0 


"B, 


J 


-k- 






,T+1 


9g + /q + eg 

E 

3=0 


j j — r 

£ £ 

r=0 fc=0 


~E t 


J 


-k- 




,fe-flGj 





which confirms ([5351) . and ([5351) . □ 
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4 Examples 



Example 4.1. Find the weighted Moore-Penrose inverse of the rational matrix 

X(s)={{s+l,s+2,s},{s,s,s+l},{s+l,s+2,s}} 

using the following weighting matrices, M\(s) and N\(s): 

Ml (s)={{s+l , s , s+l} , {s , s+2 , s} , {s+l , s , s+3}} ; 

Nl (s) ={{s+l , s+l , s+l} , {s+l , s+2 , s} , {s+l , s , s+3}} ; 

The following result is generated applying the function WPartit, implementing 
Algorithm 2.1 {see implementation details): 



WPartit [X,M1,N1] 

WEIGHTED MOORE-PENROSE INVERSE= 

/ 2s 2 (2+s) (2+s) ~ 



12+32s+33s 2 +14s 
2+5s+2s 2 



2+3s+2s 2 

2(l+a) 

(6+7s)(2+3s+2s 2 ) 2+3s+2s 2 
2s 2 (1+s) (l+s)(2+s) 
2+3s+2s 2 



s(l2+16s+5s 2 ) \ 



V (6+7s)(2+3s+2s 2 ) 



4+2s-2s^ 
12+32s+33s 2 +14s a 

s(l + s)(6 + 5s) 
(6+7s)(2+3s+2s 2 ) / 



Example 4.2. In this example we compute the weighted Moore-Penrose inverse 
of the rational matrix X(s) due to the following weights Mi(s) and Ni(s): 

X={{l/s"2,s, (s+l)/s-3},{s,s-2-l,s},{s+l,l/s,s+l}} ; 
Ml={{s+1 , s , s+l} , {s , s+2 , s} , {s+l , s , s+3}} ; 
Nl={{s+l,s+l,s+l},{s+l,s+2,s},{s+l,s,s+3}}; 

WPartit [X,M1,N1] 

WEIGHTED MOORE-PENROSE INVERSE= 

/ _ g 3 -l-s+s 5 +s 6 -1-s+s 2 



v 



2-s+s 2 + 

1 + s 
-2-s + s 2 +s 3 

_zl+nl+£^ 

-2-s + s 2 +s 3 



2+s-s 2 - 

B-S 3 +S : 

-2-s + s 2 4 



Example 4.3. If the matrices are considered in the polynomial form, then the 
function WPartPoly, implementing Algorithm 3.1, can be used to compute the 
weighted Moore-Penrose inverse of the matrix X (see implementation details): 

X={{l+s , -2+s"4 , s} , {s , -l+s , s} , {s , s , 1+s}} ; Ml=Nl={{l+s , s , s} , {s , -1+s , s} , {s , s , 1+s}} ; 



WPartPoly [X, Ml, HI] 

WEIGHTED MOORE-PENROSE INVERSE= 

1 2+2 s+s 2 -s 4 -s 5 

s — s 2 +s 5 — 1 + s + s 2 — s 3 

s T+2s 

s — s 2 +s 5 — 1 + s + s 2 — s 5 

s s (3+s — s 4 ) 

s + s 2 — s b 1 — s — s 2 +s b 



1-s-s 2 - 
-1 + 2 s + s 2 - 



Example 4.4. In this example we generate the Moore-Penrose inverse of the 
matrix X(s), known as the parameter test matrix of Hessenberg form \32\ : 

X={{s ) l,0 > 0,0},{s"2,s ) l,0,0},{s"3,s"2,s,l,0},{s"4 > s"3,s"2,s > l},{s"5,s"4 > s"3,s"2,s}} 

using identity matrices Ml(s) and Nl(s) of appropriate orders, we get: 

WPartPoly [X , Ident ityMatrix [5] , IdentityMatrix [5] ] 
WEIGHTED MOORE-PENROSE INVERSE= 

( (ifsF \ 



10 

-si 



(1 + s) 2 



V 







1 

TT+Ip 



(1+s) 
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5 Conclusion 

We extend Wang's partition method from [28] to the set of one- variable rational 
and polynomial matrices. In this way, we obtain an algorithm for symbolic 
computation of the weighted Moore- Penrose inverse of one- variable rational and 
polynomial matrices. The paper is a generalization of the paper [28] and a 
continuation of the paper [23) . Several symbolic examples are arranged. In 
partial case M = I m , N — I n we obtain the usual Moore- Penrose inverse, and 
then use test examples from [32] . Main implementation details are described as 
the appendix in the next section. 

6 Implementation details 

For the sake of completeness we describe the MATHEMATICA code which imple- 
ments Algorithm 2.1. and Algorithm 3.1. 

6.1 Rational matrix case 

Main problem in the implementation of Algorithm 2.1 is the simplification of 
algebraic expressions included. This difficulty imposes its implementation in a 
symbolic computational package. Moreover, a significant problem in the im- 
plementation of Algorithm 2.1 is the magnification of arithmetic operations. 
This problem increased by multiplicative recomputations. In view of Step 2 in 
Algorithm 2.1 , for each i £ {2, . . . ,n}, the Moore-Penrose inverse Xi(s) must 
be computed n — i + 1 times. Moreover, in view of Step 2.1 and Step 2.3, the 
pseudoinverse Xi-i(s) is needful during the computation of the values di(s) and 
bi(s). Consequently, Algorithm 2.1 requires 3(n — i + 1) recomputations of the 
Moore-Penrose inverse Xi(s), for each i € {2, . . . , n}. The total number of dif- 
ferent values that will be produced is comparatively small, but these values must 
be recomputed many times by means of relatively complicated expressions. In 
order to obviate recomputations, we use possibility of the programming pack- 
age MATHEMATICA to define functions that remember values they have found 
[311 130] . The pattern for defining a memo function is f [x_] : =f [x] =rhs. 

In order to enable simplifications of rational expressions by means of MATHE- 
MATICA function Simplify, we restrict our implementation to the set of rational 
matrices with real coefficients. 

In the beginning we describe two auxiliary procedures. 

A. The function Col[a, j] extracts j-th column of the matrix a — A(s): 

Col[a_List, j_] := Transpose [{Transpose [a] [ [j] ] }] 

B. The submatrix Aj(s) = [a\{s), ■ ■ ■ aj(s)] which contains first j < n columns 
of the matrix A(s) =A n (s) = [ai(s), • • • a„ (s)] is generated as follows: 

Adop [a_List , j_] : =Module [{m,n} , 
-Cm,n}=Dimensions [a] ; 

Return [Transpose [Drop [Transpose [a] ,-(n-j)]]] ;] 
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Step 2 of the Algorithm 2.1 is implemented in the following functions which 
remember before computed values. 
Implementation of Step 2. 1 . 

DD [a_List ,mO_List ,nO_List , i_] : =DD [a,mO , nO , i] = 
Module [{s ={}}, 
s=Simplify [A [a,mO ,n0 , i-1] .Col[a,i]] ; 
Return [s] ] 

Implementation of Step 2.2. 

CC [a_List ,mO_List ,nO_List , i_] : =CC [a.mO , nO , i] = 
Module [{s={}}, 
s=Col [a , i] -Adop [a , i- 1] . DD [a , mO , nO , i] ; 
Return [Simplify [s] ] ] 

Implementation of Step 2.3. 

B [a_List ,mO_List ,nO_List , i_] : =B [a,mO ,n0 , i] = 
Module [{nul ,ml , j ,k,nl , s={}} , 
{ml ,nl}=Dimensions [CC [a,mO ,n0 , i] ] ; 
nul=Table[0,{j,l,ml},{k,l,nl}] ; 
If [CC [a,mO,nO,i]=!=nul, 

s=Inverse [Transpose [CC [a ,m0 , nO , i] ] . mO . CC [a, mO , nO , i] ] 

. Transpose [CC [a, mO, nO, i]].mO, 
s=(Delt[a,mO,nO,i] )~(-l) . (Transpose [DD [a, mO,nO, i] ] .NK[nO,i] [[1]] 

-Transpose [NK [nO , i] [[3]]] ) . A [a,mO ,n0 , i-1] ] ; 

Return [Simplify [s] ] ] 

The following function Delt[a, mO, nO, i] computes 6{ defined in (|2.8|) . 

Delt [a_List ,mO_List ,nO_List , i_] :=Module [{s}, 

s=NK[nO,i] [[2]]+Transpose[DD[a,mO,nO,i]] .NK[nO,i] [[1]] . DD [a,mO ,n0, i] - 
(Transpose [DD[a,mO,nO,i]] .NK[nO,i] [ [3] ] +Transpose [NK [nO , i] [[3]]] . DD [a,mO ,n0 , i] ) 
-Transpose [NK [nO , i] [[3]]] . (IdentityMatrix [i-1] -A [a,mO ,n0, i-1] .Adop[a,i-l] ) 
. Inverse [NK [nO , i] [ [1] ] ] . NK [nO , i] [ [3] ] ; 

Return [Simplify [s] ] ] 

In the function NK[a, i] we find the partition (|2.2|) of the leading principal 
submatrix A^(s) of the weighted matrix N(s). 

NK [a_List , i_] : =Module [{lk , NK1 ,nkk} , 
nkk={{a[[i,i]]}}; 
If [i==l,Return[{nkk,nkk,nkk}] , 

NKl=Transpose [Take [Transpose [Take [a, i-1]] ,i-l]] ; 

lk =Transpose [{Most [Last [Take [Transpose [Take [a, i] ] , i] ] ] }] ; 
Return [{NK1 , nkk , lk}] ] ] 

Implementation of Step 1 and Step 2.4- 

A [a_List ,mO_List ,nO_List , i_] : =A [a,mO ,n0 , i] = 
Module [{b=a} , 

If[i==l, (* Compute XI (s) *) 

If [Col [a, i] ===Col [a , i] *0 , 

b=Transpose [a] [[1]] , (* al(s)=0 *) 

b=Inverse [{Transpose [a] [[i]] .mO . Col [a, i] }] . {Transpose [a] [[1]] .mO}] , (* al(s)!=0 *) 

(* Compute Xi(s), i>l *) 
b=A [a, mO, nO, i-1] -(DD [a, mO,nO,i]+ (IdentityMatrix [i-1] -A [a, mO.nO, i-1] .Adop [a, i-1] ) 

. Inverse [NK [nO , i] [[1]]] .NK[nO,i] [ [3] ] ) . B [a , mO , nO , i] ; 
b=Append[b,B[a,mO,nO,i] [[1]]]] ; 
Return [Simplify [b] ] ] 



22 



M.B. Tasic, P.S. Stanimirovic, M.D. Petkovic 



The following function starts recursive computations in Step 2: 

WPartit [a_List ,mO_List ,nO_List] : = 
Module [{m,n, i} , {m,n}=Dimensions [a] ; 
Print ["WEIGHTED MOQRE-PENRQSE INVERSE="] ; 
A[a,mO,nO,n] // MatrixForm] 

6.2 Polynomial matrix case 

We also restrict the implementation to the set of polynomial matrices with 
real coefficients. The matrix A(s) defined in (|3.1|) can be represented as the 
list {Ax, . . . ,A q +i\. The i-th column cii(s) of A(s) is the polynomial matrix 
defined in Q3.3p . and therefore can be represented by the three-dimensional list 
{di,!, . . . , a^ g+ i}, 1 < i < n. 

Col [L_List , j_] : = (* Compute j-th column from L *) 
Module [{Ll=L2={>,i}, 

For [i=l , i<=Length [L] , i++ , 

Ll=Append [LI , Transpose [L [ [i] ] ] ] ; AppendTo [L2 , Transpose [{LI [ [i , j] ] }] ] ] ; 
Return [L2] ] ; 

FrmPoly [M_List] : = C* Form the polynomial matrix of the form (3.1) *) 
Module [{L={} , i , M1=M , v , s} , 
v=Variables [M] ; 
If [v=!=0, 

s=v[[l]]; (* The matrix is not constant *) 

For[i=l, i<=Max [Exponent [M,s]] ,i++, 

AppendTo[L, Coefficient [M,s"i]] ; Ml=Ml-Coef f icient [M, s~i] *s~i] ; 
M1={M1}; 

For[i=l,i<=Length[L] ,i++, AppendTo [Ml ,L [ [i] ]]] ] ; 
If [v= !={>, Return [Simplify [Ml] ] , C* The matrix is not constant *) 
Return [Simplify [{Ml}]]] (* The matrix is constant *) ]; 

TakeFPoly [L_List , j_] : = C* Separate first j columns from L *) 
Block [{Ll={},i}, 

For[i=l,i<=Length[L] ,i++, Ll=Append [LI .Take [Transpose [L [ [i] ] ] ,j]]] ; 
Return [LI]] ; 

DopZero [L_List , i_] : = C* Complete the matrix L by zero rows *) 
Module [{Ll=L,j ,nula}, 
nula=Ll[[l]]*0; 

For [ j =1 , j <=i-Length [L] , j ++ , AppendTo [LI , nula] ] ; 
Return [LI] ] ; 

LastZeroP [L_List] : = C* Drop the last zero rows from L *) 

Module [{Ll=L,Us=True ,nul ,dl} , 
If [Ll=!={>, 

While [Us bk Ll=!={}, dl=Dimensions [LI] [ [1] ] ; 

If [Ll[[dl]]==Ll[[dl]]*0, Ll=Drop[Ll,-l] , Us=False]]]; 
Return [LI]] ; 

DDP [L_List ,M_List ,N_List , i_] : =DDP [L,M,N, i] = (* Compute d_{i,j+l} using (3.8) *) 
Module [{Y={} , gr=0 ,bb , NN , L2={> , L1=L , j , nula={}> , 
L2=ZZP[L,M,N,i-l] ; gr=Length [L] +Length [L2] ; 

nula={Table [0 , { j , 1 , gr}] } ; Ll=DopZero [LI , gr] ; NK=Col [LI , i] ; L2=DopZero [L2 , gr] ; 
For[j=0,j<gr-1, j++, 

If [(j+l)>Length[Y] ,Y= Join [Y, nula]] ; 

Y[[j+l]]=Sum[L2[[j-k+l]] .NN[[k+l]] ,{k,0,j}] ; 

] ; 

Y=LastZeroP[Y] ; Return [Y] ] ; 

CCP[L_List,M_List,N_List,i_] :=CCP[L,M,N,i]= (* Compute c_{i,j+l} using (3.9) *) 
Module [{Y=L4={} , gr=0 , NN , L1=L , L2=L3={} , j ,nula={}} , 

L2=YYP[L,M,N,i-l] ; gr=2Length [L] +Length [L2] ; nula=Table [0 , {j , 1 ,gr}] ; 
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Ll=DopZero [Ll.gr] ; NN=Col [LI , i] ; L2=DopZero [L2 ,gr] ; L4=DDP [L , M , N , i] ; 

If [L4==O,L4={0}] ; L4=DopZero [L4 ,gr] ; L3=TakeFPoly [L, ; L3=DopZero [L3,gr] ; 

For[j=0, j<gr-l,j++, 

If [(j+l)>Length[Y] ,Y=Join[Y,nula]] ; 
If [(Length [L4[[l]]] ==0) , 

Y [ [ j+1] ] =Sum [NN [ [j-k+1] ] L2 [ [k+1] ] - (L3 [ [j-k+1] ] L4 [ [k+1] ] ) [ [1] ] , {k,0 , j}] , 

Y [ [ j + 1] ] =Sum [NN [ [ j -k+1] ] L2 [ [k+1] ] -Transpose [L3 [ [j -k+1] ] ] . L4 [ [k+1] ] , {k , , j }] 

]] ; 

Return [LastZeroP [Y] ] ] ; 

VVP[L_List,M_List,N_List,i_] : =VVP [L,M,N, i] = (* Compute V_{i,j+1} using (3.11) *) 
Module [{M0=M , Ll={} , L2={} , L3={> , L4={} , L5={} , L6={} ,L7={} , 
L8={},Y={}, j ,k,r,iz,q,mq,ql,gr}, 
L2=CCP[L,M,N,i] ; L5=DDP [L , M , N , i] ; L6=NKP[N,i] [[1]] ; 
L7=NKP[N,i] [[3]] ; L8=ZZP [L ,M,N, i-1] ; 
If[L2=!={}, Ll=YYP[L,M,N,i-l] ; 

mq=Length[M0]-l; q=Length [L] -1 ; ql=Length [LI] ; 
gr=q+2*Length[Ll]+mq; Ll=DopZero [LI , gr] ; L2=DopZero [L2 , gr] ; 
M0=DopZero[M0, gr] ; iz = {}; 
For[j = 0, j < gr, j++, 
iz= Join [iz , {Sum [Sum [Sum [LI [ [j -k-r+1] ] Transpose [L2 [ [k+1] ] ] . 

M0[[r+l]],{k,0,j-r}],{r,0,j}]][[l]]}]], 

(* Else *) 

iz={};L4=Delt[L,M,N,i] [[1]] ; 

gr=Length [L4] -l+2*Length [XPP [L , M , N , i-1] ] +Length [L] -1+Length [N] -1 ; 
Ll=YYP[L,M,N,i-l] ; LI = DopZeroLLl, gr] ; L4 = DopZero[L4, gr] ; 
L5 = DopZero[L5, gr] ; L6 = DopZero[L6, gr] ; L7 = DopZero[L7, gr] ; 
L8 = DopZero[L8, gr] ; Y = {}; 
For[j = 0, j < gr, j++, 
If [Length[Dimensions[L5[[l]]]] == 1, 
Y=Join [Y , Sum [L5 [ [j -k+1] ] L6 [ [k+1] ] 

-Transpose [L7[ [j-k+1]]] LI [[k+1]] , {k,0,j}]] , 
Y= Join [Y , Sum [Transpose [L5 [ [j -k+1] ] ] . L6 [ [k+1] ] 

-Transpose[L7[[j-k+l]]]Ll[[k+l]] , {k,0,j}]]]] ; 
Y=DopZero[Y, gr] ; 
For[j=0,j<gr,j++, 

If [Length[Dimensions[L8[[l]]]] == 1, 

iz = Join[iz,{Sum[Sum[L4[[j-k-r+l,l]]Y[[k+l]] . 

{L8[[r+1]]>, {k,0,j-r}], {r , , j }] }] , 
iz = Join [iz, {Sum [Sum [L4[[j -k-r+1, l]]Y[[k+l]]. 

L8[[r+1]], {k,0,j-r}], {r ,0 , j}] }] ] ] ] ; 
Return [LastZeroP [iz] ] ] ; 

WWP[L_List,M_List,N_List,i_] : =WWP [L,M,N, i] = (* Compute W_{i,j+1} using (3.12) *) 
Module [{Y={} ,M0=M , gr=0 , iz , L0={} , Ll={} , L2={} , L3={} , L4={} , j , nula={} ,mq} , 
L2=CCP[L,M,N,i] ; L3=YYP [L,M,N, i-1] ; L0=Delt [L , M , N , i] [ [2] ] ; 
If [L2= ! ={} , iz={} ; gr=2*Length [L2] +Length [MO] -1 ; 
L2=DopZero [L2 , gr] ; M0=DopZero [MO , gr] ; 
For[j=0,j<gr,j++, 

iz=Join[iz, Sum [Sum [Transpose [L2[[j -k-r+1]]] .M0[[k+1]] .L2[[r+1]] ,{k,0,j-r}] ,{r,0, j}] [[1]]]] , 
gr=2*Length [L3] +Length [L0] -1 ; 
L4=Transpose [{Delt [L,M,N,i] [[2]]}] ; 
iz={} ; gr=Length [L4] -l+2*Length [L3] ; 
L3=DopZero [L3 , gr] ; L4=DopZero [L4 , gr] ; 
For[j=0,j<gr,j++, 

iz=Join[iz,Sum[Sum[L4[[j-k-r+l]]L3[[k+l]]L3[[r+l]] ,{k,0,j-r}] ,{r,0,j}]]]] ; 
Return [LastZeroP [iz] ] ] ; 

ZZP[L_List,M_List,N_List,i_] :=ZZP[L,M,N,i]= (* Compute Z_{i, j+1} using (3.17) and (3.6) *) 
Module [{L1=L , L2={} , L3={} , L4={} , L5={} , L6={} , L7={} , L8={} , 
M0=M,mq,rez,NNl ,q,gr ,gr2,iz, izl, j ,k,r}, 
If[i==l, (* Step 1 *) 
mq=Length [MO] -1 ; q=Length [LI] -1 ;L2={} ; 
For [j =1 , j <=Length [Col [LI , 1] ] , j ++ , 

L2=Join [L2 , Transpose [Col [LI , 1] [ [j ] ] ] ] ] ; 
If [LastZeroP [L2] ==={} , rez=L2 , 

L2=DopZero[L2,mq+q+l] ; M0=DopZero [MO ,mq+q+l] ; iz={}; 
For [ j =0 , j <q+mq+ 1 , j ++ , 
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iz=Join[iz,{Sum[Sum[L2[[j-k+l]] .M0[[k+1]] ,{k,0, j}] ] }] ] ; 
rez=LastZeroP [iz] ; If [rez=={> , rez={L2 [ [1] ] . MO [ [1] ] *0}] ] , 
OElse *) 

L4=VVP[L,M,N,i] ; iz=TET [L,M,N, i] ; L7=KSIP [L ,M,N , i] ; 
gr2=Length [iz] +2 ; iz=DopZero [iz , gr2] ; 
L4=DopZero[L4,gr2] ; L7=DopZero [L7 , gr2] ; izl={}; 
For[j=0 ) j<gr2,j++ ) 

izl=Join [izl , {Sum [L7 [ [j -k+1 , 1] ]L4 [ [k+1] ] ,{k,0, j}] }] ] ; 
rez={}; 

If [LastZeroP [iz] ==={} , iz=izl*0] ; 
For[j =0, j<gr2, 

If [Length [Dimensions [iz [ [1] ]] ] == 1, 

rez=Join [rez , {Join [{iz [ [j + 1] ] } , {izl [ [j+1] ]}]}], 
If [Dimensions [iz[[l]]] [[1]] == 1, 

rez=Join [rez, {Join [iz [ [j+1]] ,{izl [ [j+1] ] }] }] , 
rez=Join [rez , {Join [iz [ [j+1] ] , {izl [ [j+1] ] }] }] 
]]]];(* EndlF *) 
Return [rez] ] ; 

YYP [L_List ,M_List ,N_List , i_] : =YYP [L,M,N, i] = (* Compute Y_{i,j+1} using (3.18) and (3.7) *) 
Module [{Ll=L,L2=O,L3={},L4={},L5={>,M0=M,iz={},q,mq, j,k,r,gr}, 
If [i==l , mq=Length [MO] -1 ; q=Length [LI] -1 ; L2={} ; 
For [j =1 , j <=Length [Col [LI , 1] ] , j ++ , 

L2=Join [L2 , Transpose [Col [LI , 1] [ [j ] ] ] ] ] ; 
If [LastZeroP [L2] ==={} , iz=L2 , 
L2=DopZero [L2 ,mq+2*q+l] ; M0=DopZero [MO ,mq+2*q+l] ; 
L3=DopZero[Col[Ll,l] ,mq+2*q+l] ;iz={}; 
For[j=0, j<2*q+mq+l, j++, 

iz=Join[iz,Sum[Sum[Sum[L2[[j-k-r+l]] .M0[[k+1]] .L3[[r+1]] ,{k,0,j-r}] ,{r,0,j}]]]] ; 
iz=LastZeroP [iz] ] , 
(* Else *) 

L3=KSIP[L,M,N,i] ; L5=WWP [L ,M,N, i] ; 

gr=Length [L3] +Length [L5] -1 ; 

L5=DopZero [L5 , gr] ; L3=DopZero [L3 , gr] ; iz={} ; 

For[j=0,j<gr,j++, 

iz=Join [iz , Sum [L3 [ [j -k+1] ] L5 [ [k+1] ] , {k , , j }] ] ] ; 
iz=LastZeroP [iz] ] ; 
Return [iz] ] ; 

NKP [L_List , i_] :=NKP[L,i]= (* Find the partition (2.2) *) 
Module [{lk={} , NK1={} , nkk={} , Ll={} , L2={> , L3={> , L4={} , L5={}} , 
For[j=0,j<Length[L] , j++ ,nkk=Join [nkk, {{L [ [j+1] ] [[i,i]]}}]] ; 
If [i==l,Return[{nkk,nkk,nkk}] , 

For[j=0,j<Length[L] , j++ ,L1= Join [LI , {Take [L [ [j+1] ] ,i-l] }]] ; 
For [j =0 , j <Length [L] , j ++ , L2=Join [L2 , {Transpose [LI [ [j +1] ] ] }] ] ; 
For[j=0,j<Length[L] , j++ ,NKl=Join [NK1 , {Take [L2 [ [j + 1] ], i-1] >] ] ; 
L5={} ; L4={} ; L3={} ; L2={} ; Ll={} ; 

For[j=0,j<Length[L] , j++ ,Ll=Join [LI , {Take [L [ [j+1] ] ,i]}]] ; 
For [j =0 , j <Length [L] , j ++ , L2=Join [L2 , {Transpose [LI [ [j +1] ] ] }] ] ; 
For[j=0,j<Length[L] , j++ ,L3= Join [L3, {Take [L2 [ [j+1] ] ,i]}]] ; 
For[j=0,j<Length[L] , j++ ,L4= Join [L4, {Last [L3 [ [j + 1] ]]}] ] ; 
For[j=0,j<Length[L] , j++ ,L5= Join [L5, {Most [L4 [ [j + 1] ]]}] ] ; 
For [j =0 , j <Length [L] , j ++ , lk= Join [lk , {Transpose [{L5 [ [j+1] ]}]}]] ; 

]; 

Return[{NKl,nkk,lk}] ]; 

KSIP[L_List,M_List,N_List,i_] :=KSIP[L,M,N,i]= (* Compute \psi_{i,j+l} (3.21) *) 
Module [{Ll={} , L2={} , L3={} , L4={} , L5={} , iz , gr , j } , 
Ll=YYP[L,M,N,i-l] ; L2=ZZP [L,M,N, i-1] ; L3=InvNKP[N,i] ; 
L4=L3 [ [3] ] ; gr=Length [LI] +L3 [ [4] ] ; iz={} ; 
Ll=DopZero[Ll,gr] ; L4=DopZero [L4,gr] ; 
For[j=0,j<gr,j++, 

iz=Join[iz,{Sum[Ll[[j-k+l]]L4[[k+l]] ,{k,0,j}]}] ;] ; 
iz=LastZeroP [iz] ; Return [iz]] ; 

FIP[L_List,M_List,N_List,i_] :=FIP[L,M,N,i]= (* Compute \varphi_{i , j+1} (3.20) *) 
Module [{Ll={} , L2={} , L3={} , L4={} , L5={} , L6={} , iz , izl , Y , Yl , gr , gr 1 , gr2 , gr3 , j } , 
Ll=YYP[L,M,N,i-l] ; L2=ZZP [L ,M,N, i-1] ; L3=InvNKP [N,i] ; L4=L3 [[1] ] ; 
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L6=NKP[N,i] [[3]] ; L5=TakeFPoly [L, i-1] ; 
grl=Max [Length [LI] -1 , Length [L2] +Length [L] -1] ; 
LO=DopZero[FrmPoly[IdentityMatrix[i-l]] ,grl] ; 
Ll=DopZero [LI , grl] ; Yl={> ; 
For[j =0,j<grl,j++, 

Yl=Join [Yl , {Sum [LI [ [ j -k+1] ] LO [ [k+1] ] , {k , , j }] }] ] ; 
gr2=Length [L2] +Length [L] -1 ; 

L2=DopZero[L2,gr2] ;L5=DopZero[L5,gr2] ;Y = {} ; 
For[j =0, j<gr2, j++, 

If [Length[Dimensions[L2[[l]]]] == 1, 

Y= Join[Y,{Sum[{L2[[j-k+l]]}.Transpose[L5[[k+l]]] ,{k,0,j}]}] , 

Y= Join[Y,{Sum[L2[[j-k+l]] . Transpose [L5 [ [k+1] ] ] ,{k,0,j}]}] 

]]; 

gr3=L3 [ [2] ] +Length [N] ; L4=DopZero [L4 , gr3] ; 
L6=DopZero[L6,gr3] ; izl={}; 

For[j=0,j<gr3,j++,izl=Join[izl,{Sum[L4[[j-k+l]] .L6[[k+1]] ,{k,0,j}]}]] ; 
gr=grl+gr3; Y=DopZero [Y,gr] ; Yl=DopZero [Yl ,gr] ; 
izl=DopZero [izl ,gr] ; iz=0; 
For[j=0,j<gr,j++, 

iz=Join[iz,{Sum[(Yl[[j-k+l]]-Y[[j-k+l]]) .izl [[k+1]] ,{k,0,j>]>]] ; 
iz=LastZeroP [iz] ; Return [iz]] ; 

TET [L_List ,M_List ,N_List , i_] : =TET [L,M,N, i] = (* Compute \Theta_{i , j+l> (3.19) *) 
Module [{Ll={} , L2={} , L3={} , L4={} , L5={} , L6={} , iz , izl , iz2 , iz3 , gr , grl , j , kl , k2} , 
Ll=ZZP[L,M,N,i-l] ; L2=InvNKP [N, i] ;grl=L2[[4]] ; 
L2=L2[[3]]; L3=WWP [L ,M,N , i] ; L4=VVP [L,M,N, i] ; 
L5=DDP[L,M,N,i] ; L6=FIP [L,M,N, i] ; iz={}; 
gr=Length [LI] +grl+Length [L3] +Length [L4] -1 ; 

Ll=DopZero [LI ,gr] ; L2=DopZero [L2 ,gr] ; L3=DopZero [L3 ,gr] ; L4=DopZero [L4 , gr] ; 
If [L5==0 , For [j =1 , j <=i , j ++ , L5=Append [L5 , {{0}}] ] ; ] ; 

L5=DopZero [L5 , gr] ; izl={} ; iz2={} ; iz3={} ; 

For[j =0,j<gr, j++, 

izl= Join [izl , {Sum [Sum [LI [ [j -k-r+1] ] L2 [ [k+1 , 1] ] L3 [ [r+1] ] , {k , , j -r}] , {r , , j>] }] ; 
If [L5[[l]]*0==={0>, 

iz2=Join [iz2 , {Sum [Sum [ ({L5 [ [j -k-r+1] ] }L2 [ [k+1 , 1] ] ) . {L4 [ [r+1] ] } , {k , , j -r}] , <r , , J}] }] , 
iz2=Join [iz2 , {Sum [Sum [ (L5 [ [ j -k-r+1] ] L2 [ [k+1 , 1] ] ) . {L4 [ [r+1] ] } , {k , , j -r}] , {r , , j }]}]]] ; 
If [L6=={> ,L6=Table [{0} , {kl , 1} , {k2 , i-1}] ] ; 
L6=DopZero [L6 , gr] ; 
For[j=0,j<gr,j++, 

iz3=Join [iz3 , {Sum [L6 [ [j-k+1] ] . {L4 [ [k+1] ] } , {k, , j}] }] ] ; 
For [j=0,j<gr,j++, 
If [i==2 , iz=Join[iz , {{izl [ [j+1] ] }-iz2 [ [j + 1] ] -iz3 [ [j+1] ] }] , 
iz=Join [iz ,{izl [ [j + 1] ] -iz2 [ [j+1] ] -iz3 [ [j+1] ]}]]]; 
iz=LastZeroP [iz] ; If [iz=={> , iz={{0}}] ; 
Return [iz]] ; 

WPartPoly[L_List,M_List,N_List] := (* Implementation of Algorithm 3.1 *) 
Module [{mm ,nn ,k ,rez={} , Ll={} , L2={} , L3={} , Ml={} , Nl={}} , 
{mm,nn}=Dimensions [L] ; 

A=FrmPoly[L] ;Ml=FrmPoly[M] ; Nl=FrmPoly [N] ; 
For [k=l , k<=nn-l , k++ , 

Ll=ZZP[A,Ml,Nl,k] ;Print["ZZP=",Ll] ; 

L2=YYP[A,Ml,Nl,k] ; 

If [L1===L1*0, rez={Ll,{l}>, rez=SimplP [LI , L2] ] ; 
ZZP[A,Ml,Nl,k]=rez[[l]] ; YYP [A, Ml ,N1 ,k] =rez [ [2] ] ; 

L2=Sum [rez [ [2 , j]] (Variables [L] [ [1] ] ) ~ ( j -1) , { j , 1 , Length [rez [ [2] ] ] }] ; 
Ll=Sum[rez[[l, j]] (Variables [L] [[1]])~ (j-1) ,{j , 1 , Length [rez [ [1] ] ] }] ; 
Ll=VVP[A,Ml,Nl,k+l] ; L2=WWP [A ,M1 ,N1 ,k+l] ; 
rez=SimplP[Ll,L2] ; VVP [A, Ml ,M1 ,k+l] =rez [ [1] ] ; 
WWP[A,Ml,Nl,k+l]=rez[[2]] ; 

]; 

Ll=ZZP[A,Ml,Nl,nn] ; L2=YYP [A ,M1 ,N1 ,nn] ; rez=SimplP [LI ,L2] ; 
Print ["ZZP[" ,nn,"]=" ,rez[[l]]] ; 
Print ["YYP [" ,nn,"]=" ,rez[[2]]] ; 

L2=Sum [rez [ [2 , j]] (Variables [L] [[l]])~(j-l) ,{ j , 1 , Length [rez [ [2] ] ] }] ; 
Ll=Sum[rez[[l, j]] (Variables [L] [[1]])~ (j-1) ,{ j , 1 , Length [rez [ [1] ] ] }] ; 
Return [Simplify [L1/L2] //MatrixForm] ] ; 
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PolLCM [L_List] : = (* Find the least common multiple *) 
Module [{m=ml=l , j} , 

For [ j = 1 , j <=Length [L] , j ++ , 
If [Variables [L] =!={}, 

m=PolynomialLCM[m,L[[j]]] , 
m=LCM[m,L[[j]]] 

] ]; 

If [Variables [m]==0, 
If [Length [m]=!=0, 

For[j=l, j<=Length[m] ml=LCM [ml ,m [ [j] ] ] 

] ], 

If [Not [Head [m] = ! =List] , 

For[j=l, j<=Length[m] ml=PolynomialLCM [ml ,m [ [j] ] ] ], 
ml=m 

] 1; 

Return [Expand [ml] ] ]; 

SimplP[Ml_List,M2_List] := 

Module [{p , q, r , vr={} , M3=M4={} , i} , 

p=Sum[Ml [ [i+1] ] *w"i ,{i , , Length [Ml] -1}] ; 

q=Sum [M2 [ [i+1] ] *w"i , {i , , Length [M2] -1}] ; 

If [Head [q] = ! =List , r=Simplif y [p/q] , r=Simplif y [p/q [ [1] ] ] ] ; 

M3=PolLCM [Denominator [r] ] ; 

If [Variables [M3] =!={}, M4=Expand [Simplify [r*M3] ] ; 

M3=Transpose[FrmPoly[{M3}]] [[1]] ; 

M4=FrmPoly[M4] , M4=FrmPoly [r] ; M3={1}] ; 
Return [{M4,M3}]] ; 

Delt[L_List,M_List,N_List,i_] :=Delt [L,M,N, i] = (* Compute \Delta_{i , j+1} (3 . 15) , (3 . 16) *) 
Module [{Ll={} , L2={} , L3={} , L4={} , L5={} , L6={} , L7={} , gr , grO , grl , gr2 , 
rez , dell , del2 , del21 , del22 , del23 , del24 , del25} , 
Ll=YYP[L,M,N,i-l] ; L2=InvNKP [N , i] ;grl=L2[[4]] ;gr2=L2[[2]] ; 
L2=L2[[3]]; L3=NKP [N , i] [ [1] ] ; L4=NKP [N,i] [[2]] ; 
L5=NKP[N,i] [[3]] ; L6=DDP [L,M,N, i] ; L7=FIP [L,M,N, i] ; 
If [L6=={} , L6={Transpose [L5 [ [1] ] ] *0}] ; If [L7==0 , L7={L6 [ [1] ] *0}] ; 
gr=2*Length[Ll]+grl-l; Ll=DopZero [LI ,gr] ; L2=DopZero [L2 ,gr] ;dell={}; 
For [j =0 , j <gr , j ++ , dell=Join [dell , {Sum [Sum [LI [ [j -k-r+1] ] LI [ [k+1] ] L2 [ [r+1] ] , 

{k,0,j-r}],{r,0,j}]}];]; 

grO=3*Length [LI] +gr2+gr ; Ll=DopZero [LI , grO] ; 
L2=DopZero [L2 , grO] ; L3=DopZero [L3 , grO] ; L4=DopZero [L4 , grO] ; 
L5=DopZero[L5, grO] ; L6=DopZero [L6 , grO] ; L7=DopZero [L7 , grO] ; 
del2={> ; del21={} ; del22={} ; del23={} ; del24={} ; del25={} ; 
If [i==2 , For [ j =0 , j <grO , j ++ , 

del21=Join [del21 , {Sum [Sum [Sum [L4 [ [j -k-r-t+1] ] LI [ [k+1] ] LI [ [r+1] ] L2 [ [t+1] ] , 

{k,0,j-t-r}] ,{r,0, j-t}] ,{t,0,j}]}] ; 
del22=Join [del22 , Sum [Sum [Sum [Transpose [{L6 [ [j -k-r-t+1] ] }] . 

L3[[k+1]] [[l]]L6[[r+l]]L2[[t+l]] ,{k,0,j-t-r}] ,{r,0,j-t}] ,{t,0,j}]]; 
del23=Join [del23 , Sum [Sum [Sum [Transpose [{L6 [ [j -k-r-t+1] ] }] . 

L5[[k+1]] [[l]]Ll[[r+l]]L2[[t+l]] ,{k,0,j-t-r}] ,{r,0,j"t}] ,<t,0,j>]]; 
del24=Join [del24 , Sum [Sum [Sum [Transpose [L5 [ [j -k-r-t+1] ] ] . 

L6[[k+l]]Ll[[r+l]]L2[[t+l]] ,{k,0, j-t-r}] ,{r,0,j"t}] ,{t,0,j}]] ; 
del25=Join [del25 , Sum [Sum [Transpose [L5 [ [j -k-r+1] ] ] . 

L7[[k+l]]Ll[[r+l]] ,{k,0,j-r}] ,<r,0,j}]]] , 
For [ j =0 , j <grO , j ++ , 

del21=Join [del21 , {Sum [Sum [Sum [L4 [ [j-k-r-t+1] ] LI [ [k+1] ] LI [ [r+1] ]L2 [ [t+1] ] , 

{k,0, j-t-r}] ,{r,0,j-t}] ,{t,0,j}]}] ; 
del22=Join [del22 , Sum [Sum [Sum [Transpose [L6 [ [j -k-r-t+1] ] ] . L3 [ [k+1] ] . 

L6[[r+l]]L2[[t+l]] ,{k,0, j-t-r}] ,{r,0,j-t}] ,{t,0,j}]] ; 
del23=Join [del23 , Sum [Sum [Sum [Transpose [L6 [ [ j -k-r-t+1] ] ] . 

L5 [ [k+1] ] LI [ [r+1] ]L2[ [t+1] ] ,{k,0, j-t-r}] ,{r,0,j-t}] ,{t,0,j}]] ; 
del24=Join [del24 , Sum [Sum [Sum [Transpose [L5 [ [ j -k-r-t+1] ] ] . 

L6[[k+l]]Ll[[r+l]]L2[[t+l]] ,{k,0, j-t-r}] ,{r,0,j-t}] ,{t,0,j}]] ; 
del25=Join [del25 , Sum [Sum [Transpose [L5 [ [ j -k-r+1] ] ] . 

L7[[k+l]]Ll[[r+l]],{k,0,j-r}] ,{r ,0 , j}] ] ] ] ; 
del2=del21+del22-del23-del24-del25; dell=LastZeroP [dell] ; del2=LastZeroP [del2] ; 
rez=SimplP[dell,del2] ; 
Return [rez]] ; 
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(* Compute inverse of N *) 

NNinv[N_List,i_] :=NNinv[N,i]= (* Compute (3 . 33) , (3 . 34) *) 

Module [{L0={},Ll={},L2={},L3={} > L4={},L5={} > Y={},e={} > f={},g={},rez}, 
Ll=NKP[N,i] ; L2=L1 [ [2] ] ; 

If [i==l , Print ["Ninv [" ,i , "] =" , {{{1}} ,L2}] ;Return [{{{1}} ,L2}] , 
Y={}; L3=EEI[N,i] [[2]] ; L4=FFI [N , i] [ [2] ] ; L5=GII [N , i] [ [2] ] ; 
L6=EEI[N,i] [[1]] ; L7=FFI [N , i] [ [1] ] ; L8=GII [N, i] [ [1] ] ; 
gr=Length [L3] +Length [L4] +Length [LB] ; 

L3=DopZero [L3 , gr] ; L4=DopZero [L4 , gr] ; L5=DopZero [L5 , gr] ; 
For[j=0,j<gr,j++, 

Y=Join [Y , Sum [Sum [{L3 [ [ j -k-r+1] ] }L4 [ [k+1 , 1] ] LB [ [r+1 , 1] ] , {k , , j -r}] , {r , , j}] ] ] ; 
Y=LastZeroP [Y] ; e={} ; 

gr=Length [L6] +Length [L4] +Length [LB] ; L6=DopZero [L6 , gr] ; 
L4=DopZero [L4 , gr] ; LB=DopZero [LB , gr] ; 
For[j=0, j<gr, j++, 

e=Join[e,{Sum[Sum[L6[[j-k-r+l]]L4[[k+l,l]]LS[[r+l,l]] ,{k,0,j-r}] ,{r,0,j}]}]] ; 
If [LastZeroP [e] ==={} , e={e [ [1] ] } , e=LastZeroP [e] ] ; 
f=Oi gr=Length [L7] +Length [L3] +Length [LB] ; L7=DopZero [L7 ,gr] ; 
L3=DopZero [L3 , gr] ; LB=DopZero [LB , gr] ; 
For[j=0, j<gr,j++, 

f=Join[f ,{Sum[Sum[L7[[j-k-r+l]]L3[[k+l]]L6[[r+l,l]] ,{k,0,j-r}] ,{r,0,j}]}]] ; 
If [LastZeroP [f ] ==={} , f ={f [ [1] ] } , f =LastZeroP [f ] ] ; 

f =DopZero [f , Length [e] ] ; g={} ; gr=Length [L8] +Length [L3] +Length [L4] ; 
L8=DopZero[L8,gr] ;L3=DopZero [L3 ,gr] ; LB=DopZero [LB ,gr] ; 
For[j=0,j<gr,j++, 

g=Join [g , {Sum [Sum [L8 [ [ j -k-r+1] ] L3 [ [k+1] ] LB [ [r+1 , 1] ] , {k , , j -r}] ,{r , , j}] }] ] ; 
If [LastZeroP [g] ==={} , g={g [ [1] ] } , g=LastZeroP [g] ] ; 
iz=FrmPoly [FormE [e , f , g] ] ; rez=SimplP [iz , Y] ; 
Return [rez] ] ] ; 

GII[M_List,i_] :=GII[N,I]= (* Compute (3 . 27) , (3 . 28) *) 

Module [{L0=O, Ll={} ,L2=0, L3={} , L4={} , LS={} , iz, izl, iz2, gr, j ,k,r, rez}, 
LO=NNinv[N,i-l] ;Ll=NKP[N,i] ;L2=L1[[2]] ;L3=L1[[3]] ; 
L4=L0 [ [1] ] ; LB=LO [ [2] ] ; gr=Length [L4] +Length [LB] +2*Length [N] ; 
L2=DopZero[L2,gr] ;L4=DopZero [L4,gr] ; LB=DopZero [LB , gr] ; 
L3=DopZero [L3 , gr] ; izl={} ; iz2={} ; 
For[j=0,j<gr,j++, 

iz2=Join [iz2 , {Sum [L2 [ [j-k+1] ] LB [ [k+1] ] , {k,0,j}]}]] ; 
For[j=0, j<gr, j++, 

If [i<=2 , izl=Join [izl , {Sum [Sum [- (Transpose [L3 [ [ j -k-r+1] ] ] . L4 [ [k+1] ] ) 

•L3[[r+1]] ,{k,0,j-r}] ,{r,0,j}]}] , 
izl=Join [izl , Sum [Sum [-Transpose [L3 [ [ j -k-r+1] ] ] . L4 [ [k+1] ] 

.L3[[r+1]] ,{k,0,j-r}] ,{r,0,j}]] ;]] ; 

iz=iz2+izl ; 

Return [{LastZeroP [LB] , LastZeroP [iz] }] ] ; 

FFI [N_List , i_] : =FFI [N , i] = (* Compute (3.29) , (3 . 30) *) 

Module [{L0={} , Ll={} , L2={} , L3={} , L4={} , LB={} , iz , izl , gr , j , k , r , rez} , 
L0=MNinv[K,i-l] ; Ll=NKP[N,i]; L2=GII [K , i] [ [2] ] ; 
L3=L1[[3]]; L4=L0[[1]]; LB=L0[[2]]; gr=Length [L4] +Length [N] ; 
iz={}; L3=DopZero[L3,gr] ; L4=DopZero [L4 ,gr] ; 

For[j=0, j<gr, j++, iz=Join [iz ,{Sum [L4 [ [j-k+1] ] .L3[[k+1]] ,{k,0,j}]}] ;] ; 
If [LastZeroP [iz] ==={} , iz={iz [ [1] ] } , iz=LastZeroP [iz] ] ; 
Return [{-iz , LastZeroP [L2] }] ] ; 

EEI [N_List , i_] : =EEI [N , i] = (* Compute (3 . 31) , (3 . 32) *) 

Module [{L0={} , Ll={} , L2={} , L3={} ,L4={} ,L6={} , L6={} , L7={} , L8={} , 
sl={} , s2={} , iz , izl ,gr , j , k ,r ,rez} , 
L0=NNinv [K , i-1] ; L1=L0 [ [1] ] ; L2=L0 [ [2] ] ; 

LB=FFI [N , i] [ [1] ] ; L6=FFI [N , i] [ [2] ] ; L7=GII[N,i] [[1]] ; L8=GII [N, i] [ [2] ] ; 

gr=Max [Length [LI] +Length [L7] +2*Length [L8] , Length [L2] +Length [L8] +2*Length [L7] ] ; 

iz={}; Ll=DopZero[Ll,gr] ; L7=DopZero [L7 ,gr] ; L6=DopZero [L6 ,gr] ; 

LB=DopZero [LB ,gr] ; L2=DopZero [L2 , gr] ; sl={}; 

For[j=0,j<gr,j++, 

If [Length [Dimensions [LB [ [1] ] ] ] ==1 , 

sl= Join [si , {{Sum [Sum [LI [ [j -k-r+1] ] L7 [ [k+1] ] L6 [ [r+1] ] , {k , , j -r}] , {r , , j }] }}] , 
sl=Join [si , {Sum [Sum [LI [ [ j -k-r+1] ] L7 [ [k+1] ] L6 [ [r+1 , 1] ] , {k , , j -r}] , {r , , j }]}]]] ; 
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s2={}; 

For[j =0,j<gr, j++, 

If [Length [Dimensions [L5 [ [1] ] ] ] ==1 , 

s2=Join [s2 , {{Sum [Sum [L2 [ [ j -k-r+1] ] L5 [ [k+1] ] LB [ [r+1] ] , {k , , j -r}] , {r , , j }] }}] , 

s2=Join [s2 , {Sum [Sum [L2 [ [ j -k-r+1] ] LB [ [k+1] ] . Transpose [LB [ [r+1] ] ] , {k , , j -r}] , {r , , j }]}]]] ; 

iz=sl+s2; grl=Length[L7]+Length[L2]+Length[L6] ; izl={}; 

L2=DopZero [L2 , grl] ; L7=DopZero [L7 , grl] ; L6=DopZero [L6 , grl] ; 

For[j=0,j<grl,j++, 

izl=Join [izl , {Sum [Sum [L2 [ [ j -k-r+1] ] L7 [ [k+1] ] L6 [ [r+1 , 1] ] , {k , , j -r}] , {r , , j }] }] ] ; 
Return [{LastZeroP [iz] ,LastZeroP [izl] }] ] ; 

FormE [e_List ,f _List ,g_List] := 
Module [{el, f 1 ,gl , Y, Yl , Y2, i} , 
el=Sum [e [ [i+1] ] *w~i , {i , , Length [e] -1}] ; 
f l=Sum [f [ [i+1] ] *w"i , {i ,0 , Length [f ] -1}] ; 
gl=Sum [g [ [i+1] ] *w"i , {i , , Length [g] -1}] ; 
If [Head[gl]=!=List,gl={gl}] ; 
Yl=0*el; 

For[j=l, j<=Length[el] , j++, 

If [Length [f 1] ==1 , Yl [ [ j ] ] =Append [e 1 [ [ j ] ] , f 1 [ [ j ] ] ] , 
Yl[[j]]=Join[el[[j]],fl[[j]]]]]; 
If [Length [f 1] ==1 , Y2={Join [f 1 , gl] } , 

Y2={Join [Transpose [fl] [[1]] ,gl]}] ; 
Y=Join[Yl,Y2] ; 
Return [Y] ] ; 

InvNKP[L_List,i_] :=InvNKP[L,i]=Module[{LO={},rez,iz}, (* Compute (3.3S) *) 
L0=NNinv[L,i-l] ; rez=SimplP [L0 [ [1] ] ,L0[[2]]] ; 

If [i==2 , iz={{rez [ [1] ] } .Length [rez [ [1] ] ] -1 .Transpose [{rez [ [2] ] }] .Length [rez [ [2] ] ] -1} , 
iz={rez [ [1] ] , Length [rez [ [1] ] ] -1 , Transpose [{rez [ [2] ] }] , Length [rez [ [2] ] ] -1}] ; 
Return [iz] ] ; 
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